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ABSTRACT 

We study the evolution of Hernquist profile "galaxies" in the presence of self-interacting dark mat- 
ter (SIDM), where the properties of the dark matter can be parameterized by one number, &dm — 
(TdmMt/o 2 for a halo of mass Mt and break radius a. While the halos form constant density cores 
of size ~ a/2 on the core radius relaxation time scale (t rc ~ i-7td yn /o , DM) core collapse begins shortly 
thereafter and a steeper 1/r 2 central density cusp starts forming faster than predicted by 2-body re- 
laxation. The formation of the steeper central cusp is accelerated if the cooling baryons adiabatically 
compress the dark matter. The natural consequence of SIDM is to exacerbate rather than to mitigate 
astrophysical problems created by dark matter density cusps. 

Subject headings: cosmology: theory - galaxies: formation - methods:numerical - dark matter 
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1. INTRODUCTION 

A model based on the growth of small fluctuations 
through gravitational instability in a universe with cold 
dark matter (CDM) provides an excellent fit to a wide 
range of observations on large scales (3> IMpc). However 
the nature and properties of the CDM, apart from its be- 
ing cold and dark, remain mysterious. It is one of the 
major goals of cosmology to further constrain the prop- 
erties of the dark matter and determine its nature. Re- 
cently much attention has been refocused on this qu estion 
because of a suggestion by Spergel & Steinhardt (2000) 
that possible discrepancies with observations on kpc scales 
could be probes of dark matter properties, particularly 
its self-interaction rate. In its simplest incarnation, the 
self-interacting dark matter (SIDM) model provides a 1- 
parameter family of models with standard CDM as a lim- 
iting case, and it is thus interesting to examine this model 
in some detail. 

In this paper we calculate the quantitative effects of 
SIDM on the evolution of isolated halos using an N-body 
code with particle scattering. We review the astrophysical 
arguments motivating SIDM in ^ and describe our imple- 
mentation of the scattering algori thm i n ^ (our method 
is very similar to that of Burkert (2000), but uses a more 
accurate approximation to the scattering rate integral). 
The impact of self-interaction on various astrophysical pro- 
cesses is discussed in §|] and we finish with some conclu- 
sions in §[|. 

2. THE CASE FOR SIDM 

In this section we bri efly re cap the arguments which led 
Spergel & Stei nhard t ( |2000| ) to revive the SIDM mode l 
(Carlson et al. 1994 ; de Laix, Scherrer & Schaeffer 1995 ) 
since it is precisely these astrophysical observations that 
can be used to constrain the self-interaction of dark mat- 
ter. 

There is a growing consensus that CDM halos form with 
central density cusps , where p ~ r ~ 7 as r — > with 7 ~ 1 



(Navarro et al. 1997 ; Moore et al. 
Klypin et al. |1999| ; Moore et al. 



1998; Huss et al. 1999 



Kull 1999) suggest that such a core structure follows nat- 
urally from merging in collisionless systems. These re- 
sults are (in some senses) confirmed experimentally by the 
similar stellar de nsity cusps found in early-type galaxies 
(e.g. Faber et al. 1997), which also largely formed through 
the merging of collisionless matter (stars). However, such 
density cusps are inconsistent with the finite core radius 
required to repro duce the rotation c urves of dwarf galax- 



ies (Moore p94[ Flores & Primack |l99j). The problem 
appears to be limited to very low rotatio n spee d galaxies 
(v c < 10 2 kms" 1 ) as van den Bosch et al. ( |2000| ) have now 
shown that the slightly more massive LSBs (low surface 
brightness galaxies) actually require cusped dark matter 
profiles. CDM simulations also produce halos which cor- 
rectly reproduce the numbers of galaxies in clusters, but 
at first sight grossly overpredict the nu mber of satellites 
orbiting the Galaxy (Moore et al. 1999| ). 

One solution to the discrepancies is to blame the differ- 
ences on the effects of the baryons. The dwarf galaxies are 
the systems where the baryons can most easily modify the 
core struct ure of the dark matte r (e.g. Hernquist & Wein- 
berg \L992j N avarro Eke & Frenk |1996|; Gela to & Sommer- 
Larsen |1999| ; Binney, Gerhard fc Silk pOOOj ) although it is 
difficult to quantitatively estim ate the effects of feedback 
(e.g. MacLow & Ferrara 1999). Similarly, massive spiral 
galaxies are probably the best environment for the effects 
of the baryons to modify the satellite predictions of pure 
dark matt er simulations (e.g. Bullock, Kravstov & Wein- 
berg [2000 ). It is, however, unclear whether the baryons 
can provide a complete solution to the two problems. A 
second solution is to use hot rather than cold initial condi- 
tions for the colla psing dark matter (e.g. Hozumi, Burkert 
& Fujiwara pOOOl ). 

Self-interacting dark matter may be a third solution, be- 
cause it can also heat the cores of dark matter halos and 
evaporate orbiting satellites given an appropriately tuned 
interaction cross section. This idea has att racted consider- 
able atten tion ( Ostriker 2000|; Hannestad 2000]; H ogan & 
Dalcanton p000| ; Burkert poTjoj ; Firmani et al. |200C| ; Mofc 



2000a; but see Kravs- 



tov et al. [19981 ). Analytic arguments (Syer & White [1998 



Mao 2000]; Hannestad & ScEerrer 200C , Bento eTaT. 2000), 



but the detailed quantitative consequences of the theory 
have yet to be derived. We have used N-body methods, 
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discussed in the next section, to explore the formation and 
evolution of "galaxies" in a simple model. 

3. NUMERICAL METHODS 

Throughout we shall study 10 5 p article realizations of a 
Hernquist profile (Hernquist 1990) 
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which has two parameters: the total mass Mr and break 
radius a. The orbital time at the break radius a is 
t dyD = 4:n{a 3 /GM T ) 1/2 , which we shall sometimes refer 
to as the "dynamical time". With 10 5 particles we can 
probe the halo structure down to ~ 0.1a, which is suffi- 
cient for our purposes. We evolve the system using a Tree 
code described in Appendix |A|. 

In addition to the usual gravitational force, we imple- 
ment the self-interaction of the dark matter using a Monte- 
Carlo technique. First let us define some dimensionless 
units. For a galaxy of mass Mo and radius Rq (where Mo = 
Mp and Rq = a for the Hernquist model) define a density 
scale po = M /i?o, a velocity scale vo = (GAIq/Ro) 1 ^ 2 and 
a time scale to = Rq/vq. The simplest self-interaction has 
the scattering cross section per unit mass of the dark mat- 
ter, <jdm, independent of energy. In this case we obtain 
the 1-parameter family of models which we shall study 
here. The scattering rate for a particle with velocity Vq is 



T = ^ = J rf 3 vi/(vi)pCT_D M |v - Vl| 



(2) 



where /(vi) is the local distribution function of velocities. 
In dimensionless units, 



r 



dn 
di 



odm I d 3 v/(vi)/5|vo - Vi| 



(3) 



with dimensionless cross section a dm = M a DM /Ro- For 
a Hernquist model the fraction of the particles interact- 
ing per crossing time to is 0.05<tda/, the mean free path 
is a I pa dm j an d these interactions have a typical radius 
~ a/2. We shall quantify the effect of self-interaction as a 
function of a dm i an d one can then scale to cases of astro- 
physical interest by suitable choice of the basic dimensional 
parameters Mo and Rq. 

We have simulated dimensionless cross sections of odm = 
M T (Jdm/o? = 0, 0.3, 1.0, 3.0 and 10.0. For a dm = 1 
we also ran a case with 8 x 10 5 particles, doubli ng the 
spatial resolution, and a case using the Burkert (2000) 
scattering algorithm for c ompa rison. The evolution of 
Navarro, Frenk & White (1997) models should be simi- 
lar if we scale the two models to have the same central 
density distribution. For this normalization, the evolution 
of an NFW model will match that of our Hernquist models 

if a dm = 2irS c p cr i t r c aDM- 

To implement the scattering in an N-body code we need 
to discretize the calculation of the integral to compute the 
scattering rate for each particle. Call the particle whose 
rate we are calculating and particle j the particle from 
which it is scattering. We approximate the local density 
from the total mass and radius enclosed by the N = 32 
nearest neighbors which we order in increasing distance 



from particle 0. Writing the masses and velocities of the 
neighbors as rhj and Vj with j = 1, • • • , N we estimate the 
scattering rate for particle by 



T — odmPn 
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where 



Pn 



3 E J= i^ 



(4) 



(5) 



is the estimate of the local density. Note t hat this dif- 
fers from the prescription of Burkert ( [2000 ) who simply 
approximated |Av| = |vo|, thereby underestimating the 
heating rate for low velocity particles. 

For time step At the probability of an interaction during 
the time step is P n t = TAt. We adjust the time steps 
to keep the interaction probability relatively low (Pint < 
0.5), and the actual number of interactions in each time 
step is determined from the Poisson distribution with an 
expectation value of Pint- For each scattering event we 
interact with one of the N nearest particles, where the 
probability of interacting with particle j is 



Pi 



vol 



vol 



(6) 



This makes the interaction less "point-like" , but it pro- 
vides a more correct estimate of the scattering rates and 
their dependence on particle energy. We have run one 
simulation (a dm = 1) with 8 times as many particles to 
estimate the systematic error induced by this non-locality, 
and find that it is negligible (see Fig. 0a). 

Once a particle is selected for the interaction, we assume 
that the scattering is isotropic. For particles with center 
of mass velocity v c and velocity difference Av, we select a 
random direction e and assign the particles velocities of 



vo = v c + \ Av\em,j / (mo + m 3 ) 
Vj = v c - |Av|em /(m + rrij). 



(7) 
(8) 



Such elastic scattering will conserve energy and linear mo- 
mentum, but not angular momentum because of the finite 
particle separations. However, the net violation of angu- 
lar momentum conservation over all interactions will be 
zero because the particle separation vectors have random 
orientations. The finite interaction range also leads to dif- 
fusion effects, but these introduce negligible error as we 
can see by comparing to our higher resolution simulation 
(Fig- ||a). 

Finally we note that several simulations have been 
run approximating self-interacting dark matter as a fluid 
an d using numerical hydrodynamic techniques (Moore et 
al. |2000b| ; Yoshida et al. |2000|). This is a physically differ- 



ent regime from the lo w opt ical depth system proposed by 
Spergel & Steinhardt ( 2000 ). Unlike a strongly interacting 
fluid, the energy transfer in SIDM is dominated by stochas- 
tic collisions in the core which rapidly transport the energy 
into the outskirts of the halo. The fluid results are thus 
not directly applicable. The physical problem SIDM most 
closely resembles is 2-body relax ation in glob ular c lusters 
(se e Spe rgel & Steinhardt |2000|; Hannestad pOOOfc Burk- 



ert |ooc]). 
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4. RESULTS 



We considered a series of experiments comparing the 
evolution of N-body systems with and without SIDM pa- 
rameterized by a dm . The key to understanding their evo- 
lution is the relaxation time due to the scattering. It is, 
however, a different regime from normal 2-body relaxation. 
The scattering process is direct rather than diffusive and 
the structure of the system evolves on the relaxation time 
scale rather than on a large multiple of the relaxation time 
scale as in globular clusters. We hrst define a relaxation 
time for SIDM, and then compare the evolution of three 
dynamical systems with and without SIDM. We begin with 
the evolution of Hernquist profile galaxies, next we con- 
sider the effects of adiabatic compression of the dark mat- 
ter by the baryons in modifying this evolution, and finally 
we examine the collapse of a cold top-hat perturbation. 

4.1. Relaxation time 

For 2-body relaxation, the relaxation time is defined to 
be a 2 / D(vp ara ) where a is the one-dimensional velocity 
dispersion and -D(v 2 ara ) is the diffusion coefficient paral- 
lel to the particle velocity (see Binney & Tremaine 1987 ; 
§8). For SIDM, the scattering angle average of the ve- 
locity change is just l|Av| 2 so the analogy of the 2-body 
relaxation diffusion coefficient is 

D ( v2 ) = \ yp™|Av| 2 /(v )/(vi)d 3 v d 3 Vl (9) 

(10) 
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assuming that the particle velocity distributions /(v) are 
Maxwellians with one-dimensional velocity dispersions of 
a. The relaxation time is then 



1 



3a 2 



D(v 2 ) 3pa D M<7 



(11) 




Fig. 1. — (Top) The halo density profiles at t ~ t^ yn for models 
with (T dm = 0, 0.3, 1.0, 3.0 and 10.0. The dotted line shows the 
initial Hernquist profile. Note that after only one dynamical time 
the &DM — 3 and 10 cases have started to core collapse. The force 
softening is indicated by h. 
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Fig. 2. — (Top) Core radius ev&Tution asa function of dynamical 
time. The solid lines show the evolution of the core radius (tem- 
porally smoothed) for dimensionless cross sections of &dm = 0, 
0.3, 1.0, 3.0 and 10.0. Note that the standard (10 5 particles) and 
"Large-N" (8 X 10 particles) simulations for ctum = 1 show iden- 
tical evolution. (Middle) Central density evolution as a function of 
relaxation time t rc = l-7tdyn/&DM ■ As expected, the minimum 
core density occurs at t rc ~ 1. (Bottom) Precollapse. The recollapse 
of the core after reaching maximum expansion near t/t rc = 1 is 



self-similar for a time scale of t rc a 



1/2 
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where the factor of 3 arises because we considered the 
total velocity rather than a single component. The re- 
laxation time associated with the mean interaction ra- 
dius a/2, which we will call the core relaxation time, is 
t rc = 1.7idyn/o\DM- At the Hernquist break radius the 
relaxation time is 8idyn/o\DM and at the half-mass radius 
(2.4a) it is 35£dyn/o"r>M- As expected the evolution seems 
to occur on the relaxation time scale associated with the 
mean interaction radius a/2. 

Unlike globular clusters, the p cx 1/r cusped density 
distributions evolve rapidly, and core collap se occ urs after 
only 5 half-mass relaxation times (Quinlan |1996| ). The 2- 
body relaxation results will significantly overestimate the 
actual core collapse time scale because the SIDM scatter- 
ing is not diffusive. A particle undergoes one scattering 
with Av ~ v rather than slowly diffusing in energy due 
to many small scatterings. Crudely, we expect the sys- 
tem to form a constant density core of size ~ a/2 on the 
time scale t rc and core collapse faster than the 2-body re- 
laxation prediction of 5 half-mass relaxation times. After 
core collapse, the system will have a steeper p cx 1/r 2 cusp. 
This is just the behavior we see in the simulations, as we 
now describe. 

4.2. Core formation 

We first examined the formation of a core and the al- 
teration of rotation curves as a function of the dimension- 
less cross section a dm- Fig- shows the density distribu- 
tions at t ~ idyn for all five simulations. The scattering 
clearly produces a constant density core, while the simu- 
lation without any scattering roughly maintains the input 
density cusp. However, even after such a short time inter- 
val, the &dm — 10 simulation has begun its core collapse 
phase and has a higher central density than any of the 
other SIDM case s. Our results are very similar to those 
of Burkert ( [2000 ), but the higher scattering rates found 
with the better approximation to the scattering rate inte- 
gral mean that we obtain the same evolution with approx- 
imately 1 /3 the SIDM cross section. 

We non-parametrically estimate the core radius by the 
point where the density drops to 1/4 the central density. 
The central density is found by finding the region over 
which the Kolmogorov-Smirnov test provides a reasonable 
likelihood that the particle distribution is consistent with 
a constant density. This tends to provide us with an upper 
limit to the core radius, with an accuracy that is limited 
primarily by the finite particle number at very small ra- 
dius. Fig. g shows the evolution of the core radius with 
time and the central density with core relaxation time t rc . 
For comparison, the core radius due to the finite particle 
number and numerical relaxation effects starts on the force 
smoothing scale (0.05a) and slowly evolves to 0.1a, which 
represents a negligible systematic error for the much larger 
core radii created by the SIDM (recall our estimate of the 
core radius is in fact an upper limit). We also find that 
the evolution of the a dm — 1 system with 8 x 10 5 particles 
is identical to the standard simulation with 10 5 particles 
(see Fig. ||a). 

The system evolution is remarkably rapid, with a very 
narrow window between the formation of a significant core 
radius and its recollapse. The maximum core size is ap- 
proximately 0.4a and it forms on the core relaxation time 
scale t ~ t rc = 1.7tdyn/ a dm- Almost immediately af- 



ter reaching its maximum size the core begins to collapse 
again, leaving a brief window where the dark matter shows 
a large core radius. Although the formation of the core 
seems to be nearly independent of a dm when we scale the 
time by the core relaxation time scale, the subsequent core 
collapse phase is not. The rate of recollapse in units of t rc 
appears to accelerate as the cross section decreases (see 
Fig. |^b). Empirically we find the evolution of the cen- 
tral density after the reaching the maximum core radius is 

self-similar if we scale the time in units of t rc a^ M cx (J dm 
(see Fig. ||c) , which is a different scaling from that for nor- 
mal 2-body relaxation where the collapse depends only on 
the relaxation time (cx a dm)- We hypothesize that the 
difference arises because the scattering is not a diffusive 
process. In each scattering Aw ~ v and particles escape 
the core in a single scattering, leading to faster evolution. 

4.3. Adiabatic Contraction 

Very careful tuning of the SIDM cross section and ages 
might allow dwarf galaxies to have softened cores today 
without undergoing core collapse. However, we must ex- 
amine the effects of the baryons on the dark matter pro- 
file to estimate the appropriate cross section, as the cool- 
ing baryons will adiabatically compress the SIDM dark 
matter and significantly modify the evolution of the sys- 
tem. The adiabatic compression both raises the scattering 
rates and reduces the relaxation time scales in three ways. 
First, by raising the densities and reducing the dynami- 
cal time scale, the relaxation time drops. Very crudely, if 
we compress the dark matter by a factor of /, the scat- 
tering rate jumps by a factor of / 7 ^ 2 , producing a cor- 
responding reduction in the relaxation time. A second, 
more subtle factor is the difference in how density cusps 
and density distributions with cores are adiabatically com- 
pressed. Standard adiabatic compression models for a 1/r 
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Fig. 3. — The combined effects of SIDM and adiabatic compres- 
sion. The SIDM cross section is a = (solid) or 1 (dashed) and the 
adiahatically compressed curves include the baryonic contribution 
(Eq. n3). Curves of increasing central density are just before the 
compression (f = idyn)i just after the compression (t = 3td yn ) and 
some time after compression (t = 4.0t c i yn )- The force smoothing 
scale is 0.01a. 
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density cusp will maintain a 1/r density cusp aft er the 
compression (see Quinlan, Hernquist & Sigurdsson 1995 ) 
while the compression of a density distribution with a finite 
core radius produces a steeper (1/r 3 / 2 to 1/ r 9 / 4 ) c usp de- 



pending on the distribution function (Young 1980). Thus, 
if the SIDM produces a core radius before the adiabatic 
compression phase, the adiabatic compression produces a 
steeper than normal central density cusp. In 2-body relax- 
ation, steep c usps e volve more rapidly than shallow cusps, 
and Quinlan ( 1996| ) found that where a 1/r density cusp 
reaches core collapse in 5 half-mass relaxation times, a 
1/r 2 density cusp reaches core collapse in 1 half-mass re- 
laxation time. Hence the steeper central cusp created by 
adiabatic compression of the SIDM core further acceler- 
ates the ultimate evolution of the system to a central 1/r 2 
density cusp. Third, the existence of the baryonic matter 
provides a component of the gravitational potential whose 
evolution is not driven by scattering. 

We investigate this with a single, somewhat contrived 
experiment. We put 20% of the mass of the profile in 
an analytic Hernquist potential with a corresponding re- 
duction in the masses of the dark matter particles. We 
modeled the adiabatic compression by running the system 
for £dyn, and then slowly changing the external potential 
to 

v r + a i 

between i<jyn and 2id yn , finally fixing the external poten- 
tial for the further evolution. We used Mj, = 0.2Mt, = 
0.1a, and reduced the force smoothing to 0.01a. The final 
external potential was chosen to have the same circu lar ve - 
locity profile as a Kuzmin disk (Binney & Tremaine 1987), 
and underestimates the amount of compression co mpare d 
to the standard approximation (Blumenthal et al. 1986|) . 
It corresponds to a baryonic density 



p(r) 



3M b 



47r (r 2 + a 2 ,) 5 / 2 ' 



(13) 



At the end of the compression phase, the central dark mat- 
ter density has risen by a factor of ~ 10. We ran the 
compression sequence with either a dm = or 1. 

Fig. |^ shows the dark matter density and rotation curve 
profiles of the adiabatically compressed systems both at 
the end of the compression phase (3idyn) and the end of 
the experiment (4idyn)- At the start of the compression, 
the scattering has begun to establish the typical SIDM core 
structure. The adiabatic compression then produces sig- 
nificantly different central density profiles. The no scatter- 
ing simulation maintains a roughly p cx 1/r profile shifted 
upwards in density, while the a dm = 1 simulation has 
formed a steeper, nearly p tx 1/r 2 profile with a higher 
central dark matter density than if there had been no scat- 
tering, just as we expected from the differences in how the 
two profiles would be adiabatically compressed. In the 
subsequent evolution, the no scattering simulation main- 
tains a stable profile while the a dm = 1 simulation has a 
rising central density. The system appears to proceed di- 
rectly into the core-collapse phase of its evolution despite 
the fixed baryonic potential. 

4.4. Spherical Collapse 



Our previous experiments assumed that no significant 
scattering occurred during the initial collapse of the dark 
matter halo so that the initial conditions contained a 
p cx 1/r density cusp. We also considered how the for- 
mation of the halo is altered in the SIDM model by con- 
sidering the very simple case of the collapse of a top-hat 
overdensity with vacuum boundary conditions both with 
and without the DM scattering. These simulations are 
meant only to illustrate the primary effects of DM scatter- 
ing on collapsing perturbations, and not as realistic models 
of the formation of galaxies or clusters in a cosmological 
context. 

We generated a 10 5 particle uniform density sphere with 
velocity dispersion 1% of the virial velocity and allowed 
it to collapse. Without particle scattering the collapse 
formed a profile akin to the Hernquist profiles of the previ- 
ous sections, with some of the particles being ejected from 
the system. By fitting a Hernquist profile to the result- 
ing system we identified the scattering rate appropriate to 
a dm — 1 and re-ran the simulation from the same initial 
conditions including scattering. 

The simulations indicate that the final equilibrium pro- 
file is similar with and without scattering, however at fixed 
dynamical time the scattering run always has a higher cen- 
tral density. We postulate that the scattering provides a 
mechanism, similar to dissipation, which allows the system 
to attain higher phase space densities than are possible in 
a purely collisionless/dissipationless collapse. 

5. CONCLUSIONS 

We combined an N-body code with a Monte Carlo model 
of scattering to examine the effects of self-interacting dark 
matter (SIDM) on the properties of dark matter halos. 
Our models form a one parameter family with "standard" 
CDM as a limiting case. 

We first examined the time evolution of Hernquist mod- 
els with 1/r central density cusps in the presence of SIDM . 
Our simulations are similar to those of Burkert (2000), 



but use a better approximation to the scattering inte- 
grals. The models form a constant density core with a 
radius ~ 0.4a on the core radius relaxation time scale 
t rc = 1.7tdyn/<5 dm, where the SIDM cross section per 
unit mass is odm — o~dmMt / a 2 ■ The initial growth of 
the core radius is very rapid because the relaxatio n tim e 
is much shorter in the inner regions (see Quinlan 1996). 
Shortly after reaching this maximum core size the core 
begins to shrink. We find that the recollapse proceeds 
on a faster time scale than expected from 2-body relax- 
ation. The central density evolution evolves on time scales 



1 /2 i i 

of trc^DM ^ ® dm rather than t rc oc <Jd M , probably be- 
cause SIDM is not a diffusive process when the mean free 
paths are large. Thus, the dark matter briefly has a large, 
finite core radius which could be fine tuned to address the 
problem of dwarf galaxy rotation curves at the expense of 
predicting core collapsed dark matter halos with p cx 1/r 2 
cusps in most other galaxies. In particular, it is proba- 
bly impossible to use SIDM to eliminate the dark matter 
cusps in the dwarf galaxies while preserving them in the 
low surface brightness galaxies. 

Even with SIDM, however, we must also consider the 
role of the baryons in modifying the structure of the dark 
matter halo. In particular, as the baryons cool and form 



-1/2 
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a disk they ad iabat ically com press the dark matter (Blu- 
menthal et al. 198E; Dubinski 1994). The adiabatic com- 
pression raises the central dark matter density, increases 
the SIDM scattering rate, and reduces the time scale for 
core collapse. Moreover, adiabatic compression of a den- 
sity profile with a finite c ore rad ius produces a steep cen- 
tral density cusp (Young 198C) which relaxes to form a 
p oc l/r 2 density cusp even faster than t he sha llower 
p oc l/r cusp we considered initially (Quinlan 1996| ). As a 
result, the development of a core radius due to the SIDM 
scattering is first reversed by the compression, and then 
core collapse begins on a far shorter time scale. Such a 
further acceleration of the time scale for producing a final 
density distribution with a steep p oc l/r 2 density cusp 
will make it even harder for SIDM to maintain a core ra- 
dius in the dark matter profile for periods comparable to 
the presumed lifetimes of dwarf galaxies. 

All of the above simulations supposed that the halo col- 
lapsed and virialized before scattering became important, 
which is appropriate to "low" scattering rates. Finally, 
we considered the effects of SIDM on the collapse of cold 
top-hat perturbations with a large enough that scattering 
and virialization could occur at the same time. We found 
that the scattering allowed higher phase space densities at 
fixed dynamical time than could be found in the collision- 
less systems. 

It thus appears as though SIDM exacerbates rather than 
solves the "central density cusp problem" . Turning this 
around we can say that astrophysics almost certainly re- 
quires SIDM cross sections small enough to avoid signifi- 
cant scattering over the age of the universe for all galaxies, 
or mean free paths > 1 Mpc. 
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APPENDIX 

TREE N-BODY CODE 

We have used a new implementation of a Tree code to 
evolve our dark matter halo. Since this code is new we 
briefly describe it here, and also explain how we generated 
the initial conditions. 

The basic algorithm is described in Barnes & Hut ( 1986] ). 
Space is partitioned into an oct-tree structure by bisection 
in the three spatial dimensions. Starting at the root node, 
which is the entire volume of the simulation, the tree de- 
scends until each node contains at most one particle. The 
force on any given particle from all of the other particles 
is computed by a tree walk. The speed-up over direct 
summation is obtained by imposing an opening criterion 
during the tree walk. For cells sufficiently far from the 
particle in question, the entire mass distribution within 
that cell is approximated by a point at the center of mass 
of the cell and with mass equal to the sum of the particle 
masses within the daughters of that cell. 

Though the code does not require it, we have used equal 
mass particles throughout. The size of the computational 



box grows so as to always encompass all of the parti- 
cles. Rather than a Plummcr potenti al we use a spline 
softe ned fo rce (Monaghan & Lattanzio |1985| ; Hernquist & 
Katz |989| Springel, Yoshida & White |2000[). For most of 



the runs we use h — 0.05a (the force was therefore exactly 
l/r 2 beyond h), but for the adiabatic compression runs 
we lower this to 0.01a. Very roughly this corresponds to 
a Plummer law smoothing e ~ h/3, although a Plummer 
law gives 1% force accuracy only beyond lOe. 

With 10 5 particles it is finite particle numbers which 
limit our resolution, not force accuracy. For profiles with 
almost any central cusp, the relax ation time is less than 
lOidyn for radii r < 0.05a (Quinlan 1996). Our numerical 
results with differing particle numbers and force accuracies 
support the scalings in Fig. 2 of Quinlan ( 199(f ). Only for 
the adiabatic compression runs, where more particles are 
concentrated at small radius and there is a stable external 
potential, do we gain by lowering h. 

The time integration is done with a second order leap- 
frog method. The time step is chosen to be the smaller 
of O-Sy^/i/aniax, where a max is the maximum acceleration 
on the particle, or 2h/v where v is the particle velocity. 
We have shown that this time step accurately integrates 
binary orbits, and in combination with the opening crite- 
rion described below conserves total energy to much better 
than a percent over 5 dynamical times, except during the 
core collapse phase where the error rises to a few percent. 

Finally for the tree walk we do not simply use the geo- 
metrical criterion advocated by Barnes & Hut. Rather we 
have augmented it with a mod ification of the method de- 
scribed in Springel et al. ( [2000 ). Starting with the known 
accelerations from the last time step we descend the tree. 
We open an internal node if the partial acceleration from 
the quadrupole moment of that node exceeds a times the 
old acceleration, where a is a tolerance parameter. Specif- 
ically if I is the size of a cell containing mass M whose 
center of mass is a distance r from the particle in ques- 
tion, we open the node if 



Ml 2 > a|a i d |r 4 



(Al) 



This ensures that the relative force accuracy is 0(a) while 
the resu lting tree walk is significantly faster (Springel et 
al. 2000 ). The tree walk is also used to define the neigh- 
bor list for the density and scattering calculations. To 
ensure that all nearby cells are opened so the neighbor list 
is "complete" we augment the opening criterion to addi- 
tionally require that i < Or for cells closer than a. This 
provides a good estimate of the density except in the lowest 
density regions where the scattering probability is anyway 
small. Throughout we have used a = 2% and = 0.4, 
which for a Hernquist profile results in a 90th percentile 
relative force error of 5 x 10~ 3 , i.e. 90% of the particles 
have force errors less than this. 

We set up our initial conditions as a Hernquist ( 1990 ) 
profile 

, . Mbaio a 1 

Ph(t) = — A2) 

2ir r (r + ay 

where a is the scale-length. First we choose the position 
of a particle with a radial probability distribution propor- 
tional to r 2 p(r) and an isotropic f . We impose an upper 
radial distance of r max = 50a at which point the den- 
sity is ~ 10~ 5 of the density at r = a. Given the posi- 



tion we then calculate the magnitude of the velocity from 
the known distribution function f(E) for this model us- 
ing an acceptance-rejection method and again distribute 
v isotropically. This procedure is iterated N p /2 — 5 x 10 4 
times. Each particle and a reflected counterpart at — f 
with velocity — v is added to the list to obtain a realiza- 
tion of the halo with N p particles. The use of 'reflected' 
initial conditions ensures that the center of mass position 
and velocity of the halo are both initially zero. 
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